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We use a Fundamental-Measure density functional for hard board-like polydisperse particles, 
in the restricted-orientation approximation, to explain the phase behaviour of platelet colloidal 
suspensions studied in recent experiments. In particular, we focus our attention on the behavior of 
the total packing fraction of the mixture, n, in the region of two-phase isotropic-nematic coexistence 
as a function of mean aspect ratio, polydispersity and fraction of total volume 7 occupied by the 
nematic phase. In our model, platelets are polydisperse in the square section, of side length a, but 
have constant thickness L (and aspect ratio k = L/{a) < 1, with (a) the mean side length). Good 
agreement between our theory and recent experiments is obtained by mapping the real system onto 
an effective one, with excluded volume interactions but with thicker particles (due to the presence 
of long-ranged repulsive interactions between platelets). The effect of polydispersity in both shape 
and particle size has been taken into account by using a size distribution function with an effective 
mean-square deviation that depends on both polydispersities. We also show that the bimodality of 
the size distribution function is required to correctly describe the huge two-phase coexistence gap 
and the nonlinearity of the function 7(77), two important features that these colloidal suspensions 
exhibit. 

PACS numbers: 61.30.Cz, 61.30.Hn, 61.20.Gy 



I. INTRODUCTION 

Colloidal suspensions of mineral or viral anisotropic 
particles interacting via short-ranged repulsive forces ex- 
hibit a phase behavior with entropy-driven phase transi- 
tions between their liquid-crystalline phases. The nature 
of these phases strongly depends on particle geometry. 
The rod geometry in mineral or viral particles favours 
the formation of uniform phases, i.e. isotropic (I) and 
nematic (N) phases [l| , and also of the layered smectic 
(S) phase [2-4]. In the case of the plate geometry, the I 
and N phases are usually stabilized at low particle con- 
centration Q , and the N phase requires a relatively high 
aspect ratio (and thus particle anisotropy) Q. In addi- 
tion, as the volume fraction increases, there may appear 
a transition to the columnar (C) phase 0,0]. 

The plate geometry is more versatile as regards the 
type of liquid crystalline phases it may induce @. Re- 
cently it has been shown that colloidal suspensions of 
some plate-like, mineral charged particles, can also stabi- 
lize the smectic phase 0, [13 ■ The colloidal particles are 
usually polydisperse in their sizes (diameter and thick- 
ness) and shapes (rod or plate geometry or different par- 
ticle cross-sections), and it was found that polydisper- 
sity causes phase behavior in these systems to be much 
more complex due to phenomena such as size se greg ation, 
fractionation and multiple phase coexistence For 
example, polydisperse rod-plate mixtures exhibit demix- 
ing between I and different N phases (with the former 



populated by particles with less anisotropy), and also 
up to four coexisting phases (some of them nonuniform) 
may exist at high densities No trace of biaxial N 

phases was found in these mixtures. However, recent ex- 
periments on board-like colloidal particles could find this 
elusive phase [l2|. Also I-N coexistence with density in- 
version has been observed [ill, [TiJ (with the I phase being 
the densest phase). 

The theoretical modeling of these kind of mixtures 
turns out to be a difficult task. Density functional theory 
(DFT), which is based on the local density distribution, 
has been successful in the description of bulk and inter- 
facial phase behavior of hard spheres and other fluids of 
anisotropic particles [l5j . but it becomes hard to imple- 
ment for the case of polydisperse mixtures. This is due 
to the huge number of degrees of freedom on which the 
local density now depends: not only on spatial and (for 
anisotropic particles) orientational coordinates, but also 
on a number of polydispersity variables. Despite the in- 
creased difficulty, some theoretical calculations on poly- 
disperse mixtures of freely rotating hard rods in the On- 
sager approximation have been performed PHH?]]- These 
calculations confirm the phase behavior found in experi- 
ments as regards the broadening of the I-N coexisting gap 
and the size- fractionation phenomenon. In this respect, 
it would be interesting to extend the recently proposed 
Fundamental-Measure DFT for freely rotating hard disks 
[l8| to the calculation of phase behavior in polydisperse 
platelets. 
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Monte Carlo simulations on polydisperse infinitely thin 
hard-platelet fluids have also been carried out [l^|. These 
results again show the dramatic effect of polydispersity 
on phase behavior in hard-platelet fluids (mainly segre- 
gation driven by particle size), as compared to that in 
one-component fluids, already simulated in the 90s for 
the cut-sphere geometry [2(| . In the latter work the bulk 
phase diagram was traced. For aspect ratios k = L/a < 1 
the N and, at higher densities the C, phases were found 
to be stable (here L is the thickness and a the particle 
diameter). In the range 0.1 < k < 0.2 the C phase gives 
rise to the I and eventually for k ~ 0.2 to the cubatic 
phase, and finally for k > 0.3 only the I and the solid 
phase are stable. Similar phase behavior was found in 
Monte Carlo simulations of hard oblate spherocylinders 
(2l| . where two different crystals are stable (tilted for 
k < 0.45 and aligned for k > 0.45); the cubatic phase is 
always unstable. The practical difficulties in implement- 
ing DFT calculations associated with polydispersity can 
be circumvented by considering discrete particle orienta- 
tions, as in the Zwanzig model where the main axes of 
the particles, one for uniaxial and two for biaxial geome- 
tries, point along one of the three Cartesian axes [22j ■ In 
the framework of this model the phase diagrams of poly- 
disperse hard rods [23j and rod-plate mixtures [24| have 
been calculated. 

Within the same approach, the effect of polydisper- 
sity on the stability of the biaxial N phase of hard 
board-like biaxial particles [25[ has recently been stud- 
ied. The phase diagram of the one-component limit of 
this fluid was recently obtained [26] using a DFT based 
on Fundamental-Measure theory for hard biaxial paral- 
lelepipeds [22J. Finally, the same model has been applied 
to the study of interfacial properties of binary mixtures 
confined by external potentials (28| . 

Recently a systematic experimental study of the phase 
behaviour of polydisperse platelets in suspension has 
been presented |30|. Particles were synthesized by hy- 
drothermal methods, and further exfoliation (through 
chemical treatment with TBA molecules the pristine 
Zirconium-Phosphate (ZrP) crystals are delaminated 
into single layers) [3(| ■ The novel feature of this fluid 
is that polydispersity is in the platelet diameter, with a 
strictly constant thickness. Also the shape of the particle 
cross sections is polydisperse, with most particles hav- 
ing hexagonal geometry. The work was focused on the 
study of the I-N transition for different polydispersities 
and mean aspect ratios. The aim of the present article is 
to theoretically understand the results presented in Ref. 
[30l | . The experimental results and the conclusions we 
obtain from our theoretical model can be summarized as 
follows: 

(i) Highly polydisperse platelet mixtures exhibit a huge 
I-N coexistence gap which cannot be explained by simply 
assuming a unimodal distribution function for particle 
diameters. However, if one considers a bimodal distribu- 
tion function, the coexistence gap can be explained via 
a demixing mechanism. We remark that bimodality may 



not be apparent through a direct visual inspection of the 
bimodal distribution function. Note that the effect of bi- 
modality in platelet thickness (not diameter) on phase 
behavior has already been studied in Ref. [1J|, but in 
that case the bimodality in the size distribution function 
is clearly seen from the size histogram, which is not the 
case in Ref. [30l |. 

(ii) The repulsive character of colloidal platelet interac- 
tions in the experimental system is incorporated through 
an effective platelet thickness L c fi which is much higher 
than the thickness of the real platelets. The effective 
thickness is chosen to guarantee a reasonable description 
of experimental data by the theoretical model. 

(iii) In order to account for the shape polydispersity 
and to adequately describe the experimental findings, the 
polydispersity coefficient (mean square deviation of the 
size distribution function) should be taken much higher 
than that given in Ref. [30( . In the present study we have 
used a DFT for hard board-like polydisperse particles 
with square cross sections of width a and constant thick- 
ness L, in the oblate particle regime k = L/a < 1. The 
DFT is based on the FMT for the hard-parallelepipeds 
in the restricted-orientation approximation (27| . 

The article is organized as follows. In Sec. |TT] we 
present the theoretical model, making special emphasis 
on the implementation of the I-N coexistence calcula- 
tions (Sec. Ill A[l . the size distribution function used to 
model polydispersity (Sec. Ill B[) and the effect of particle 
shape polydispersity on the effective size polydispersity 
of the mixture (Sec. Ill Cp . Sec. IIIII presents the results, 
and is divided in two sections: Sec. IIII Al which is de- 
voted to the phase behavior of the mixture assuming a 
unimodal size distribution for diameters, and Sec. IIII Bl 
which presents results obtained with bimodal distribu- 
tions. Finally we draw some conclusions in Sec. IIVI 



II. THEORY 

The theory we use is based on a density functional for 
hard board-like particles, formulated in the restricted ori- 
entation approximation (the so-called Zwanzig model). 
Particles are polydisperse in the side length a of the 
square section, but their thickness L is fixed. The mean 
aspect ratio k = L/(o~) (with (a) the mean side length) 
is less than unity, k < 1, so that we are in the oblate- 
particle regime. The main quantities that describe our 
model are the set of density distribution functions p v (a) , 
where v = {x, y, z} refers to particles with their main 
axis parallel to the v Cartesian axis; each of the three 
orientations can be considered to correspond to a dif- 
ferent species, and the fluid can therefore be treated as 
a three-species mixture. In the following sections we de- 
scribe the theoretical formalism we have used to calculate 
the isotropic (I)-nematic (N) coexistence in polydisperse 
mixtures that possess unimodal/bimodal size distribu- 
tions. 
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A. Coexistence calculations 

If a fraction 7 of the total volume of the system V is 
occupied by a nematic (N) phase in coexistence with the 
isotropic (I) phase, then the density distribution func- 
tions p v>s {a) of the two phases, s = {I,N} should fulfill 
the lever rule (conservation of the total number of parti- 
cles): 



Po(v) =J^2p»,n(<7) + (1-7)^/^,1(0-). 



(1) 



functional minimization of @ with respect to p 1/tS (a), 
together with the definition (HJ), provide the moments at 
equilibrium [23[: 



p$ = po I da 



where 



h{a)a a e' 



7 j2 e-^N V) + 3(1 _ l)e -^r ] (°) 1 



(5) 



The density distributions of all species are the same for 
the I phase, p x .j(a) = p y ,i(a) = p z ,i{a) = pi(<r), while 
for the uniaxial N phase we have p x ,N(a) = Py,isi{a) = 
P_l,n(c) and p z .^{a) = P|| ,n (c) ; here we take the ne- 
matic director to be parallel to the z-axis. The density 
distribution function po(a) = poh{a) {parent distribution 
function) is a product of the mean parent number density 
Po and the size distribution function h{a) which fulfills 
the normalization condition J dah(a) — 1. Note that 
h(a) has units of [Length] -1 . As po is the number den- 
sity having units of [Length] -3 the distribution functions 
Po(a) and p v . s (a) have units of [Lenght] -4 . 

For the whole system I+N (with < 7 < 1), we de- 
fine a Lagrange functional from the free-energy density 
in reduced thermal units, <E> = /3F/V (with F the free 
energy, j3 = 1/kT, k Boltzmann constant and T the tem- 
perature), as 

$ [{p*,.}] = 7* [{p*,n}] + (1 - 7)* [{Pu,i}] 



dap (a) 



Po(<r) -1^2p v ,n((t) 



(2) 



where po(a) are the Lagrange multipliers that guarantee 
the constraints ([T]). Note that po(cr) is just the scaled 
with (3 chemical potential of species v of width a in each 
one of the coexisting phases, i.e. f3p, VtS (a) = po(a), V 
v = x,y, z and s = I, N. <E> is split in ideal 



$id [{P^s}} 



J da ^2 Pu,s (a) [In p v<a (a) 



I] 



(3) 



and excess 4> exc parts. In our treatment, we obtain the 
latter from Fundamental-Measure Theory (FMT). In the 
FMT formalism $ oxc ({p" s }) is a function of a finite num- 
ber of moments of the distribution function. The latter 
are defined as 



pZ 



J dap VtS {a)a a , ,0 = 0,1,2. 



(4) 



The expression for the function <3> exc ({/3" s }) [see Eq. 
(|A5[) ] is obtained in the Appendix from the scaled particle 
theory (SPT), the uniform limit of the fundamental mea- 
sure free-energy density functional (27| . The constrained 



fl,,(cxc)/ \ \ " d^cxc a 
a=0 Op v .s 



(6) 



are the excess part of the chemical potential of species v 
of width a in the phase s. Note that this is a quadratic 
polynomial in a whose coefficients are in turn functions 



of the moments p^) [2 



We use the notation p^ = p^j 



and p[°* c '(a) = p\Ji°' '(a) W. Thus the set of nine 
equations (|5|) (which guarantee the equality of chemi- 
cal potentials of all species in both phases) are solved 
self-consistently for the moments p^i, while the other 
quantity to be determined, po, is found by imposing the 
condition of mechanical equilibrium, i.e. the equality of 
pressures [see Eq. (|A8[) ] in both phases: 



(cxc) 



Pi 



Po,{p{ Q) })=P N (po,{^}) 



(7) 



The fluid pressure can be found from (|A8[) . The cloud- 
I-shadow-N coexistence, corresponding to a situation 
where the system volume is occupied by the I phase ex- 
cept for a coexisting, infinitesimal amount of the N phase, 
is obtained by taking 7 = in the expressions above. The 
case of the shadow-I-cloud-N coexistence, which corre- 
sponds to the opposite case (i.e. the N phase occupying 
the whole system volume but in coexistence with an in- 
finitesimal amount of the I phase) , is obtained by taking 
7=1. 



B. Length distribution function 

In the present study we choose a probability distribu- 
tion function h{a) which is, in the general case, bimodal: 



x (a 
h(a) = — h — 
a\ \ax 



a 2 \<Ji 



(8) 



Here x is the molar fraction when the fluid is strictly a 
binary mixture; otherwise x can be regarded as a parame- 
ter which controls the relative heights of the two maxima, 
located at a\ and 02 • The function ho(u) is selected to 
be 



ho(u) = Cu v e 



v„-Au p 



(9) 
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where the constants C and A are calculated from the 
normalization conditions J °° duho(u) — J °° duuho(u) = 



1. Thus we find 



A 1 ^ = 



T[(v + 2)/p] 



c = 



(10) 



with r(x) the Gamma function. For fixed p the param- 
eter v controls the polydispersity, while p controls the 
decay of the distribution at large a (note that for p = 1 
a Schultz distribution is obtained). All these definitions 
guarantee the normalization J °° dah(a) = 1. For the 
first moment we find J °° daah(cr) = x<j\ + (1 — x)oi. 
The one-component limit is recovered by setting x = 1 
and consequently (a) = a\. Defining the polydisperse 
coefficient for ho(u) as the mean-square deviation, Ao = 
yj (u 2 )o/(u)q — 1, where (u a )o = J °° duu a ho(u), we find 
that the polydispersity coefficient f or the bimoda l dis- 
tribution function h(a), i.e. A = J (a 2 ) / (a) 2 — 1 [with 
(a a ) = J °° dacr a h(a)) results in 



A = ^\J<J 2 A 2 + x(l - x)(a 2 - en) 2 , 



(11) 



where we have defined a a = x<r" + (1 — x)a^ ■ As a 2 > a 2 
for x 7^ or 1, we find A > Ao. For a fixed parameter r = 
<J<zjo\ > 1, the polydispersity coefficient as a function of 



x, A(x), has a maximum at x* 



with value 



A(x*) 



2^ 



1 



(12) 



For example, setting Ao = 0.3, and for r = 1.5, 2, 2.5 
and 3, we find A(x*) = 0.368, 0.476, 0.579 and 0.673 re- 
spectively; we then see that the bimodality dramatically 
increases the effective polydispersity A of the mixture. 

To present the results in the following sections we use 
the packing or volume fraction 77,5(7) with s = I, N, which 
is a function of 7, and is defined through the zeroth mo- 
ment of the distribution function as 773(7) = Ps \y)o\L. 
Specifically we will use for the presentation of results 
the values 7/1 = 371(0) = po(0)a 2 L and t?n = t?n(1) = 
po(l)(j 2 L, i.e. the cloud-I and cloud-N packing fractions. 
Also we use the total packing fraction of the polydisperse 

mixture 7/(7) = 777*1(7) + (1 - 7) r /i(7) = A)(7) cr i% tne 
latter equality being a consequence of the lever rule. Fi- 
nally we will use the length distribution functions cor- 
responding to coexisting phases defined as h^ s '(a) = 
Po 1 Si/ Pv, s { a )i ( s = I>N). Again, using the lever rule, 
we have ^\a) + (1 - l)h^(a) = h{a). 



types of polydispersities constitutes a highly nontrivial 
task, the usual procedure is to map the real particles 
onto effective particles of fixed shape but with an effective 
polydispersity in their sizes. In the following we describe 
how the effective polydispersity can be calculated in our 
particular system in which the main quantity that gov- 
erns phase behavior is the particle area of the transverse 
section (we note again that the thickness L is constant). 

We consider a system made of a collection of parti- 
cles of fixed thickness L and different cross sections. To 
be more precise, we suppose the latter to have the form 
of regular polygons inscribed in circumferences of dif- 
ferent diameters 2R and also with different number of 
edge-lengths n. Thus our system is polydisperse in the 
variables R and n, the former controlling the size poly- 
dispersity, while the latter controls the particle shape. 
The cross-sectional area of these particles can be calcu- 
lated as A n {R) = — sin ( — J . Suppose the polydis- 



2 V n , 

perse coefficient [mean square deviation of the probability 
distribution function h(R)\ is A, and we define the prob- 
ability to find a polygon with n edge- lengths as p n . Now 
we map our system onto an effective one, monodisperse 
in the number no of edge-lengths, but with an effective 
polydispersity A e g. We define the mapping as 



where 

(A n ) h = I dRh(R)A n (R) 



n) hi 



n(R) 2 (l + A 2 ) 



(13) 



2vr 



(14) 



We take (R)h = (R)h o! i, i-e. the mean radii that fol- 
low from the distribution functions h(R) and h c g(R) are 
exactly the same. Defining the coefficient 



<1 = 



I]»Pn77sin(27r/n) 
no sin(27r/no) 



(15) 



we find that the effective polydispersity coefficient can be 
found as 



A eff = v/g(l + A2)-l. 



(16) 



We have used the following expression for the probability 

Pn- 



P n = (j ^^P n -^-Pr\ n>4, (17) 



C. Polydispersity in particle shape 

In most colloidal suspensions of anisotropic particles, 
polydispersity exists not only in particle size but also in 
particle shape. As the inclusion in the theory of both 



with the triangular shape excluded, which is usually the 

case in experiments. These probabilities fulfill the nor- 

00 

malization condition p n = 1 . Once we fix the mean 

n=4 

number of edge- lengths (n) , the value of p, a function of 
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FIG. 1: Effective poiydispersity A e fi as a function of the mean 
number of edge-lengths (n) for no — 6, Ao = 0.3 and k = 9. 
Inset: probability p n as a function of n for k — 9 and (n) = 8. 



FIG. 2: The cloud-I and cloud-N coexisting packing fractions 
771 and 77N as a function of poiydispersity Ao for values of re 
equal to: (a) 0.2, (b) 0.1, (c) 0.02 and (d) 0.01. 



(n) and fc, has the form p = — — - -. The polydis- 

(ti) k o 

persity in the number of edge-lengths can be quantified 
through the coefficient 

^iw-- l =i—^—v < i8 » 

Thus the number k controls the poiydispersity. In Fig. 
[T]we plot the effective poiydispersity coefficient A e g as a 
function of (n) for the case uq = 6, Ao = 0.3 and k = 9. 
As can be seen from the figure, A c g is a monotonically 
increasing function of (n) and it can reach values above 
0.5. A particular example of p n (for k = 9 and (n) = 8) 
is plotted in the inset. 

III. RESULTS 

We have carried out coexistence calculations following 
the procedure described in Sec. Ill Al First we consider 
a polydisperse mixture with unimodal distribution func- 
tion (x = 1) and Gaussian tail (q = 2), and vary the poly- 
disperse coefficient in the range < Ao < 0.75. The re- 
sults are plotted in Fig. [2]in the r/ s ~Ao phase diagram for 
four different values of the aspect ratio: k = Ljo\ = 0.2, 
0.1, 0.02, and 0.01. 

The main conclusions we can draw from these results 
can be summarized as follows: (i) the two-phase region 
is broadened as poiydispersity is increased (in agreement 
with other theoretical and experimental results) ; (ii) this 
effect is enhanced as re is lowered, and (iii) the behaviour 
of 77N as a function of Ao changes with re: for relatively 
high values of re it is a decreasing function, whereas for 
low enough re it becomes an increasing function for large 
Ao- 

In a recent experiment [30(, suspensions of polydis- 
perse platelets were prepared from exfoliation of pris- 



tine a-ZrP crystals using TBA molecules. The result- 
ing platelets were found to be polydisperse in diameter 
and also in shape (although most particles have approx- 
imately hexagonal geometry), with a constant thickness 
equal to 26.8 A. Different aqueous suspensions were pre- 
pared with particles of mean aspect ratio [as measured 
by dynamic light scattering (DLS)] ranging from 0.001 to 
0.01. The poiydispersity coefficients of the suspensions 
were estimated using Dynamic Light Scattering (DLS), 
and all of them were found to be in the range 18%— 39%, 
with most samples having a value of about 30%. Samples 
were divided into three different sets: A, B and C. Sam- 
ples corresponding to set B were obtained from nematic 
fractionation of an original suspension followed by dilu- 
tion and, consequently, this set has the lowest poiydis- 
persity coefficients. The other two sets, A and C, result 
from the original synthesis and exfoliation of the pris- 
tine crystals. All platelets have negative surface charges 
which are partially neutralized by the positive charges of 
the TBA molecules, thus creating effective dipoles. Non- 
neutralized charges and dipoles cause the pair-interaction 
between two platelets to be long-ranged and repulsive. 



A. Unimodal length distribution 

With the aim of modeling the fluid, we mapped a col- 
lection of repulsive, shape- and diameter-polydisperse 
platelets onto effective polydisperse board-like particles 
interacting through excluded volume. To properly take 
into account the effect of long-ranged repulsive interac- 
tions, the effective thickness of particles has to be larger 
than the thickness of the actual colloidal platelets, and 

the effective aspect ratios n c g = = re— — are ob- 

(cr) L 

tained from the real re by scaling by a factor / = — ^— > 1. 
As poiydispersity increases, / should also increase due to 



FIG. 3: The cloud coexisting packing fractions r^N as a func- 
tion of the aspect ratio K e ff for polydispersities fixed to (a) 
Ao = 0.3, and (b) 0.5. Open circles, filled circles and o pen 
squares correspond to the experimental results in Ref. [30] 
for sets A, B and C (adequately rescaled with factors / = 5, 
3 and 9, respectively). 



the presence of platelets with large surface area which, as 
discussed above, should repel each other more strongly. 

This effect is in fact obtained with the model, as shown 
in Fig. [21 where the cloud packing fractions yyi N are 
plotted as a function of K c g for two values of polydisper- 
sity, Ao = 0.3 in (a) and Ao = 0.5 in (b), and using a 
unimodal length distribution h(a) with q = 2. Also in- 
cluded in the figures are the experimental results from 
[30j with / set to 5, 3 and 9 for samples A, B and 
C, respectively. These values were chosen to ensure a 
proper agreement between theory and experiment (note 
that no least-square optimization was attempted) . Sam- 
ples in set B are less polydisperse, and consequently / is 
smaller. As can be seen from the figure, these samples 
are relatively well described by our model using Ao = 0.3 
(as in the actual samples), except for those experimen- 
tal points with the two higher values of K e ff. However, 
for samples in set A, better agreement is obtained with 
Ao = 0.5. Note that polydispersity in the experimen- 
tal samples is in diameter and also in shape which, as 
discussed in Sec. Ill C[ demands that the effective poly- 
dispersity of a single-shaped model be higher. Finally, 
samples in set C, which are characterized by huge coex- 
istence gaps, are not well described by unimodal length 
distributions. The multimodality of the distribution is 
probably behind this behaviour (see following section). 

Fig. [4jshows the percentage of total volume 7 occupied 
by the N phase as a function of the total packing fraction 
r\ of the mixture for those values of K e s corresponding to 
the sets A [Fig. Ufa)] and B [Fig. [3(b)]. In the former 
case solid lines are results from calculations with a uni- 
modal distribution using q = 2 and Ao = 0.5, while in 
the latter case Ao = 0.3. It is clear that when Ao = 0.3 
the function 7(77) is practically a linear function [see Fig. 
[4] (b)], while for Ao = 0.5 some nonlinearity is already 
apparent, a trend which is more pronounced in experi- 




FIG. 4: Percentage of total volume 7 occupied by N phase as 
a function of total packing fraction n for values of K e ff corre- 
sponding to the experimental sets A, panel (a), and B, panel 
(b). Rhombi, squares, circles, triangles and stars are used to 
show the experimental results from (30| as K e s is increased. 
Solid lines correspond to the theoretical results obtained from 
a unimodal length distribution with q = 2 and (a) Ao = 0.5, 
(b) Ao = 0.3. Theoretical results using a bimodal distribu- 
tion are represented by dashed curves. Values of parameters 
that better describe experimental data represented by stars 
and triangles are as follows. Panel (a): Ao = 0.5, q = 2, 
(Ti/iej = 10, (T2/i e H = 17.5 and x = 0.7 for stars, and 
eri/Z/eff = 20, <72/£eff = 26 and x = 0.45 for triangles. Panel 
(b): Ao = 0.3, q = 2, ai/L eS = 20, a 2 /L cH = 32 and x = 0.78 
for stars, and ai/L e g = 20, U2/L e a = 32.6 and x = 0.65 for 
triangles. 



ment. Coexistence gaps from theory and experiment are 
similar, except for the highest value of K G g- in set A [stars 
in Fig. [3J (a)], and for the two higher values of K e g in set 
B [triangles and stars in (b)]. A possible reason for these 
deviations will also be explained in the following section. 



B. Bimodal distribution 

In this section we demonstrate that the multimodality 
of the length distribution function can explain both the 
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existence of a huge I-N coexistence gap, and the strong 
nonlinearity of 7(77). These are two of the main features 
present in the experiments of [3(| for those samples with 
higher values of K e g, i.e. samples in sets A and B. For 
set C we show below that the bimodality is crucial to 
adequately describe experimental results. The origin of 
this behavior is the coupling between the fractionation 
effect, typical of polydisperse mixtures, and the demix- 
ing phenomenon that occurs in multicomponent mixtures 
of particles with sufficiently different lengths. In binary 
mixtures with very asymmetric species, entropy forces 
the system to segregate into two phases of very different 
composition, and consequently the coexistence density 
gap is very large compared to that in one component 
systems or in polydisperse mixtures with unimodal size 
distributions. 

Dashed lines in Fig. @] are the theoretical predictions 
for 7(77) in the case of the experimental sets A and B 
[triangles and stars in panels (a) and (b), respectively]. 
Calculations were based on a bimodal parent distribution 
function h(cr), as described in Sec. Ill BL and the corre- 
sponding functions for the different cases are shown in 
Figs. [5] (a) and (c) (solid lines). Values for the parame- 
ters in h(a) were chosen so as to optimise agreement with 
the experimental data [stars in Figs. U] (a) and (b)]. As 
can be seen from Figs. the parent distribution func- 
tion seems to be unimodal in both cases, A and B, even 
though compositions close to 50% were chosen (x = 0.70 
and x = 0.78, respectively). The high polydispersity 
(Ao = 0.5 and 0.3, respectively) creates a large overlap 
region between the two distribution functions centered at 
(Ti and a-2 [see Eqn. ©], which results in the absence of 
a second maximum near ai . 

The combined effect of fractionation and demixing can 
also be observed in Fig. [5] from the shape of the distri- 
bution functions for the shadow-N phase [dashed curves 
in panels (a) and (c)]. Note that the fraction of platelets 
with lengths a ~ <j\ is much lower than that in the co- 
existing (parent) I phase, but the fraction for lengths 
a > a 2 increases dramatically. As a result, there appears 
a shoulder in the distribution function, and its decay for 
large a is much slower. The opposite effect occurs for the 
shadow-I distribution functions [solid lines in panels (b) 
and (d)[: now the I phase is rich in platelets with a ~ 01 , 
while the distribution function decays much faster and 
the polydisperse mixture is poor in large platelets. 

As mentioned above, the most dramatic disagreement 
between the theoretical calculations based on unimodal 
distributions and the experimental results corresponds to 
samples in set C (see Fig. [3]). One possible reason for 
this disagreement is the decay rate of the parent distribu- 
tion function, controlled by the parameter q. In order to 
check this, we have implemented our calculations using 
a truncated unimodal distribution with q = 1 (Schultz 
distribution), Ao = 0.5 and the same values of K e g as in 
Fig. [21 The results, represented by means of dashed lines 
in Fig. [(2a), give a broad coexistence gap (note that for 
vanishingly small values of 7, the total packing fraction 




a a 



FIG. 5: Length distribution functions ft' I,N ' (cr*), as a function 
of reduced length a* = cr/ai, of coexisting I (solid curves) and 
N (dashed curves) phases. Panels (a) and (c) correspond to 
cloud- I-shadow-N coexistence (i.e. 7 = 0), while (b) and (d) 
refer to cloud- N-shadow-I coexistence (7 = 1). The parent 
distribution function, which coincides with that of the cloud- 
I (for 7 = 0) or cloud-N (for 7=1) cases, is bimodal with the 
following parameters: panels (a) and (b), Ao = 0.5, q = 2, 
(Ti/Leff = 10, <72/£eff = 17.5 and x — 0.70; panels (c) and (d), 
A = 0.3, q = 2, ai/L cH = 20, a 2 /L cS = 32 and x = 0.78. 

rj rapidly increases from 771, which is shown by stars), but 
they fail to reproduce the experiments. Also, the nonlin- 
earities of the curves 7(77) are not correctly described. 

In the same figure, the theoretical results using bi- 
modal distribution functions with q — 2 and Ao = 0.5 are 
also plotted. The parameters <7j/_L e ff an d x were chosen 
to guarantee a reasonable agreement between theory and 
experiment so that now both the coexistence gaps and 
also the nonlinear behavior of 7(77) are reproduced. Fig. 
[SJb) shows the bimodal parent and shadow-I, N coexist- 
ing length distributions (corresponding to those param- 
eters which better describe the experimental points of 
Fig. Ufa) represented by circles). Note that the bimodal 
parent distribution function looks unimodal. Also, the 
shadow-N distribution exhibits a plateau in the range of 
scaled lengths 7-12, while the shadow-I is highly localized 
about the value 1. 

To illustrate how a strong nonlinearity in 7(77) can 
emerge from the coupling between polydispersity and bi- 
modality, Fig. [7J a ) shows 7 as obtained from calcula- 
tions using a bimodal distribution with q = 2, Ao = 0.3, 
x = 0.85, <J\/L c s = 20 and different values of 0-2/ L e g. As 
the latter is increased, the nonlinearity becomes stronger 
up to a point where there appears a well-defined loop in 
7(77) [see inset in Fig. EJa)]. We should mention here 
that the presence of this loop is not related with a triple 
I-N1-N2 or I1-I2-N coexistence, a fact we have proved by 
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FIG. 6: (a) Percentage of total volume occupied by N phase 
as a function of total packing fraction r\ for values of K e g cor- 
responding to the experimental samples in set C (30l |. Rhom- 
bus, squares and circles are used to represent the experimental 
values for increasing K e g- Dashed lines: results for the same 
value of K e ff but using a truncated Schultz distribution func- 
tion (q = 1) with Ao = 0.5. Stars indicate values of 771 and rjn 
for each K e ff. Solid lines: results from a bimodal distribution 
functions with q — 2 and Ao = 0.5; the other parameters 
specifying the degree of bimodality we give only for the curve 
that better describes the experimental data shown with cir- 
cles. They are: ai/L e g — 4, (T2/ 'L e g = 12 and x = 0.9. (b) 
Length distribution functions for the latter set of parameters. 
Solid, dashed and dotted lines correspond to parent, shadow- 
N and shadow-I distribution functions, respectively. Zoom of 
(b) shows second maximum in shadow-N distribution func- 
tion. 



solving the set of equations for the moments and pres- 
sures of the different coexisting phases and checking that 
they always converge to solutions corresponding to FN 
coexistence. Fig. [7J (b) shows the distribution functions 
h^' N \a) for the three different FN coexistences occur- 
ring for 7 = 0.2175 [symbols in inset of panel (a)] inside 
the loop. Clearly, as 77 is increased, the function h^ N '(cr) 
becomes less peaked at o\ but more peaked about 02 and 




FIG. 7: (a) 7 versus 77 in logarithmic scale as obtained from 
a bimodal distribution with q = 2, Ao = 0.3, x = 0.85, 
ai/L e g = 20 and <J2/L e g — 80 (dotted curve), 90 (dashed 
curve) and 100 (solid curve). Zoom in (a) is a detail of the lat- 
ter case showing the loop in the function 7(77). (b) Three pairs 
of coexisting distribution functions h^' N \a*) corresponding 
to the points shown in (a) for 7 = 0.2175. Solid and dashed 
lines represent I and N phases, respectively, while different 
colors represent different values of r\ (black, red and green in 
the order of increasing rj). 



with a slower decay for large a, while the function h^(a) 
remains practically the same. The presence of this loop 
is clearly related with the density inversion phenomenon. 
As can be seen in Fig. [7J (b), the moments of the I and 
N coexisting distribution functions fulfill the inequalities 

> Pi^ an d < Pi^ [due to the defect (excess) of 
platelets of width close to <j\ (o^) in the coexisting N 
phase with respect to the I phase]. As we have defined 
the total packing fraction 77 as being proportional to the 
zeroth moment, and this moment is lower for the N phase, 
the density inversion is produced. 

To better understand this behavior, we resort to the 
lever rule ([T]). Dividing the whole equation by po and 
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integrating over a, we find 

1 = 7 (po)^n(po) + [1 - 7(A))] Ax(po), 



(19) 



where we have defined A s = J dah^ (a) as the total area 
under the curve h^(cr) (s = I,N). Both these areas and 
7 are functions of the parent number density po- From 
(UHJ) we find that the derivative of 7 with respect to po is 
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A'^A! - I) + A[(l - A N ) 
(A N - AjY 



(20) 



In those polydisperse fluids where the parent distribution 
function is strictly unimodal, density inversion does not 
occur and we have that A\ < 1 and A-^ > 1. Taking 
the latter inequalities into account and the fact that po 
goes from p$ to p Q N ^ (the values of the parent number 
densities corresponding to the cloud-I-shadow-N coex- 
istence, with Acioud-i = 1, and the cloud- N-shadow I 
coexistence, with A c i ou d-N = 1, respectively), we have 
most likely that A[ < and A^ < 0. Thus we find from 

Eq. (HOI) that 7' > in the whole interval \p§\ p Q N) ]. 
When the parent distribution function is bimodal, den- 
sity inversion could occur. For the latter situation we 
have the opposite scenario: A\ > 1 and A-^ < 1. Thus 
we could conclude that A[ > and A^ > and then we 
obtain again 7' > [see Eq. ((20)) ] in the whole range of 
Pq. However for strong bimodality [when the two peaks 
of h(a) are well visible as shown in Fig. (b)], the sign 
of A' N could change from positive to negative giving rise, 
for certain values of po, to 7' < 0. To elucidate the condi- 
tions necessary for having a negative sign of 7', we resort 
again to Eq. (|2"0"|) . From the latter it is easy to show 
that, if 7' < 0, we obtain the condition 



— In|l- A N I > — \n\l-Aj\ 
dpo dpo 



(21) 



In Fig. [8] we plot the functions S(po) = — — In 1 1 — A s \ 

dp 

when h{a) is unimodal (a) and bimodal (b) (that corre- 
sponding to the results shown in Fig. [7]). We can see 
that while the condition ([2~Tj) is not fulfilled for any po 
for the unimodal h(a), there exists, with a bimodal h(a), 
a range of po (shaded in the figure with grey color) for 
which this condition is fulfilled. It is easy to show, us- 
ing Eq. HU), that the condition (|2"Tj) is equivalent to the 
following inequality 



714,1 > (1-7)4- 



(22) 



Thus we conclude that 7' < when the rate of change in 
the area under the curve ftW (a) weighted with the factor 
7 is greater than the corresponding rate of change of area 
of h^ T '(a) weighted with 1 — 7. We can see from Fig. [7J 
that the distribution functions h' N ' (<r) corresponding to 
the three values of po shown (p^ < p^ < p 3 ^ ) are dra- 
matically different. These values are well inside the range 
of po where A' N < [see Fig. [5] (b)[. While the first peak 
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FIG. 8: 



d 
dpo 



In 1 1 — j4i,n| (dashed and solid lines correspond 



to I and N phases, respectively) as a function of r\ = poa\L 
for h(a) taken as unimodal (a) and bimodal (b), the latter 
corresponding to the results shown in Fig. [7J The shaded 
region in (b) shows the interval of r\ where the condition (|2ip 
is fulfilled. 



centered at a ~ a 1 decreases, the second peak centered 
at a ~ o"2 increases. But these changes do not compen- 
sate each other, resulting in the net lowering of the total 
area under the curve as po increases. We can also see 
that (a) practically remains the same for these three 
values of po- A small increment of the first peak with po 
makes A[ > 0. Finally, for po = p^ the condition (f2"2")l 
is fulfilled and we have at this point that 7' < 0. 

To end this section, we show how the behaviour of the 
cloud coexisting packing fractions 771^ as a function of 
« c ff = L c ff I o\ changes with the presence of bimodality. 
These packing fractions are shown in Fig. E](a) as a func- 
tion of K e fi; results from unimodal and bimodal distribu- 
tion functions are represented by dashed and solid curves, 
respectively. A measure of how this function behaves is 
dln?7LN 



the parameter r = 



Note that if 771^ 



T cir 



for 



d In K e e 

small K ff, we obtain r ~ a, i.e. r is a measure of the 
local power-law dependence of 771. n as a function of K e g. 
In the limit K c g ~ 0, it can be shown that a ~ 1 for 
unimodal size distributions and for both I and N curves, 
with the I curve having t > 1 and the N curve r < 1 (the 
latter deviating much more from unity) [see Fig. Hfb)]. 
As polydispersity increases, this behavior is reached for 
lower values of K e ff- The parameter r as a function of 
n e f[ is plotted in panel (b). From this figure we can see 
that the bimodality dramatically increases the value of 
r corresponding to the cloud-I coexistence curve, while 
for the cloud-N coexistence curve the effect is the oppo- 
site for small enough «, but to a lesser extent. Therefore, 
some caution should be exercised in extracting the power- 
law dependence of 771^ with K e g by simply measuring the 
slope for small values of the aspect ratio. 
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FIG. 9: (a) Cloud I and N packing fractions j/^n as a function 
of K e g = L e g/<Ji for unimodal (dashed curve) and bimodal 
(solid curve) distribution functions. Parameter values are set 
as q — 2 and Aq = 0.5 for both distributions, while for the 
bimodal distribution x = 0.85 and oijo\ = a + bn c g, with a 
and b chosen such that cr^/cri = 1 for n a g = 0.2 and 02 /&i = 3 
for K e ff — 0.01. (b) r parameter, calculated from 771, n shown 
in (a) (see the text for definition), as a function of K e ff- 



IV. CONCLUSIONS 



ameter polydispersity (obviously platelets with large sur- 
faces -their number depending on the width and the tail 
of the length-distribution function-, and consequently 
with more charges and dipoles, create much more re- 
pulsive effective interactions). Using this procedure, the 
three experimental samples, separated into three distinct 
sets A, B, and C (the major difference between them be- 
ing their polydispersity) were mapped onto effective hard 
platelets with three different effective thicknesses. The 
hard model used was based on the restricted-orientation 
approximation for polydisperse board-like particles and 
a Fundamental-Measure density functional was used. 

Our theoretical results agree with experiment for sets 
A and B, except for samples with the highest aspect ra- 
tios. For set B (samples obtained via fractionation) we 
used a polydispersity coefficient Ao = 0.3, approximately 
the same value measured in experiments, and the the- 
ory correctly describes the experimental phase behaviour. 
For samples in set A, agreement is reached for Ao ~ 0.5, 
which is higher than the experimental value. However, 
we have noted (Sec. Ill C| that, since colloidal particles 
are also polydisperse in shape, once we choose the par- 
ticle geometry the value of the effective polydispersity 
should be higher. This is in fact the second important 
result of our work. 

The last result is related to the theoretical modeling of 
samples in set C. We have shown that, in this case, a uni- 
modal size distribution cannot describe the experimental 
phase behavior correctly as regards the huge density gap 
of the I-N coexistence and the strong nonlinearity in the 
percentage of volume occupied by the N phase as a func- 
tion of total volume fraction. However, a bimodal size 
distribution adequately describes the gap and the non- 
linearity; note that bimodality may not be apparent by 
direct visual inspection of the distribution function. 

In summary, both fractionation and demixing phenom- 
ena are important to explain the experimental results. 
These conclusions are expected to remain valid even if a 
more exact theoretical treatment, including free particle 
orientation, could be implemented. 



We conclude by highlighting three results of the 
present work. The first concerns the repulsive charac- 
ter of pair-interactions between the colloidal platelets 
of Ref. j3(j. The specific interactions between nega- 
tively charged platelets partially neutralized with TBA 
molecules (thus creating surface dipoles) in an aqueous 
solvent is difficult to model. In particular, hydration 
layers mediate platelet interactions. The experimental 
platelet volume fraction for the I-N transition is much 
higher than that predicted by excluded-volume-based 
models for given mean aspect ratio, implying the presence 
of highly repulsive effective interactions between parti- 
cles. It is then reasonable to map particles onto effective 
hard-core platelets with a larger effective thickness so as 
to obtain comparable results between theory and experi- 
ment. The effective thickness mainly depends on the di- 



Appendix A: Free energy following the SPT 

In this section we derive, following the SPT formal- 
ism, the expression for the free-energy density of a poly- 
disperse Zwanzig fluid made of hard board-like parti- 
cles. The fluid consists of a collection of board-like par- 
ticles with square polydisperse cross-section a and con- 
stant thickness L. The main particle axes point along 
one of the three Cartesian axes x,y,z. The microscopic 
variables describing the fluid are the density distribu- 
tion functions with v labeling the different species 
(u = {x,y,z}). 

The work to insert a scaled particle of dimensions Ai<7 
(scaled with the width parameter Ai) and X2L (scaled 
with the thickness parameter A2) and pointing along the 
v direction in the polydisperse Zwanzig fluid can be cal- 
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culated as 



PW U {X U X 2 ) = - In 



1 - 



^2 J da'p s (a') 



xK (cxcl) (Aia, X 2 L,a',L)] 



(Al) 



where the excluded volume between the particles v and s 
(the respective labels of particle orientations) , the latter 
with dimensions a' x a' x L, has the expression 

V^ cX \Xxa,X 2 L,a',L) = 11 [a UT (X 1 a,X 2 L) 

+a ST (a',L)], (A2) 

where we have defined a UT (a, L) = a + (L — a)5 UT , with 
8 UT the Kronecker delta. 

Following the SPT, the excess part of the chemical po- 
tential of the species v can be calculated as a sum of two 
terms. The first one is the second-order Taylor expan- 
sion of W V (X\, A2) around the point (0, 0) and evaluated 
at (1,1). The second one is the product of the fluid pres- 
sure P and the particle volume v = La 2 (i.e. the ther- 
modynamic work to open a cavity of dimension v) . Thus 
we have 



p^(a) = Wv (0,0)+J2 



cW„ (0,0) 



d 2 W„(0,0) 



Pv. 



(A3) 



2 ^ dXidXj 

Taking into account the thermodynamic relations 

/SJ~ /V 
dap, u {a)p v {a) - T/V, p„(a) = — 

(A4) 

with J-~[{p„}], p v and V respectively the free-energy 
density functional, the chemical potential of species v 
(which is splitted in the ideal and excess part: p v [u) = 
p( ld )(er) + p,^ K °\a)), and the system volume, we obtain 



the expression for the excess part of the free energy den- 
sity $ cxc = /3T CXC /V: 



$cxc = -"o ln(l - n 3 ) + 



ni • n 2 n 2x n2yn 2z . 
1 - n 3 (1 - n 3 ) 2 



which coincides with the uniform limit of the excess part 
of the free-energy density for a general inhomogeneous 
fluid following the FMT [13. The weighting densities 

n a are functions of the moments p^' of the distribution 
function p v {<j): 



Pi a) = 



J dap v {a)a a , a = 0,1,2. 



(A6) 



Their expressions are 



(A7) 



Finally, the pressure can be calculated from (|A5|) as 
<9<E>exc _ n ni • n 2 2n 2x n 2y n 2z 



pP 



dn 3 1 - n 3 (1 - n 3 ) 2 (1 - n 3 ) 3 



(A8) 
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